This document demonstrates the analyses reported in Gilmore, Thomas, and Fesi (2015).
In a typical set-up, we have the following file directory structure:
project-name/R # home for the .RProj, .R, .Rmd, and .html files. project-name/data # .csv files project-name/figs # .pdf or *.png files
# File paths
dir.figs <- "../figs"
dir.data <- "../data"
fn.data <- "moco-3-pattern-child.csv"
fn.egi <- "egi.csv"
# knitr options to save figures to dir.figs and create 12x6 figures at 300 DPI
knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path=paste(dir.figs, "/", sep=""))
Source libraries
# Libraries
library(ggplot2)
library(dplyr)
library(png)
library(gridExtra)
# Significance threshold
p.thresh = 0.0005
Load MOCO dataframe.
df.moco.sm <- read.csv(paste(dir.data, fn.data, sep="/"))
Analysis proceeds by analyzing the results of the low-order integer harmonics of the global motion coherence modulation (on/off at 1.2 Hz) separately then the 1F2 or dot update rate (24 Hz).
Load EGI channel positions.
df.egi <- read.csv(paste(dir.data, fn.egi, sep="/"))
# Load topo ears & nose from file
img <- readPNG(paste(dir.figs, "topoplot.png", sep="/"))
# Add a raster
topo_ears_nose <- rasterGrob(img, interpolate=TRUE)
# separate frame for annotations
featured <- c(62,75,106,107)
df.annotate <- data.frame(chan=featured, xpos=df.egi$xpos[featured], ypos=df.egi$ypos[featured], Harmonic=c("1F2", "1F2", "1F1", "1F1"))
pl.egi <- ggplot( data=df.egi, aes(x=xpos, y=ypos, label=chan)) +
geom_text(size=4) +
coord_fixed() + # 1:1 aspect ratio
scale_y_continuous("", breaks=NULL) + # omit y axis
scale_x_continuous("", breaks=NULL) + # omit x axis
geom_point(data=df.annotate, aes(x=xpos, y=ypos, shape=Harmonic), size=12) +
scale_shape(solid = FALSE) + # open shapes
annotation_custom(topo_ears_nose, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
theme( panel.background = element_rect(fill=NA),
legend.position = "none") # blank background and no legend
pl.egi
To estimate main effects for age.
this.harm = "1F1"
df <- df.moco.sm %>%
filter(Harm == this.harm)
maov.omnibus <- manova(data=df, formula = cbind(Sr,Si) ~ AgeDays + Pattern*Speed*as.factor(Channel))
summary(maov.omnibus)
## Df Pillai approx F num Df den Df
## AgeDays 1 0.000003 0.0564 2 32254
## Pattern 2 0.000006 0.0517 4 64510
## Speed 2 0.000024 0.1920 4 64510
## as.factor(Channel) 127 0.140261 19.1548 254 64510
## Pattern:Speed 4 0.000031 0.1236 8 64510
## Pattern:as.factor(Channel) 254 0.038468 2.4904 508 64510
## Speed:as.factor(Channel) 254 0.019932 1.2783 508 64510
## Pattern:Speed:as.factor(Channel) 508 0.022329 0.7169 1016 64510
## Residuals 32255
## Pr(>F)
## AgeDays 0.9451
## Pattern 0.9950
## Speed 0.9427
## as.factor(Channel) < 2.2e-16 ***
## Pattern:Speed 0.9983
## Pattern:as.factor(Channel) < 2.2e-16 ***
## Speed:as.factor(Channel) 2.234e-05 ***
## Pattern:Speed:as.factor(Channel) 1.0000
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
This analysis revealed non-significant effects of AgeDays, Pattern, and Speed, but a significant effect of Channel, \(F(254, 32255)=19.1548, p<2.2e-16\) qualified by significant Pattern by Channel, \(F(508, 32255)=2.4904, p<2.2e-16\), and Speed by Channel, \(F(508,32255)=1.2783,p<2.2e-5\) interactions.
this.harm = "1F2"
df <- df.moco.sm %>%
filter(Harm == this.harm)
maov.omnibus <- manova(data=df, formula = cbind(Sr,Si) ~ AgeDays + Pattern*Speed*as.factor(Channel))
summary(maov.omnibus)
## Df Pillai approx F num Df den Df
## AgeDays 1 0.000009 0.1336 2 31102
## Pattern 2 0.000001 0.0106 4 62206
## Speed 2 0.000005 0.0403 4 62206
## as.factor(Channel) 127 0.122587 15.9913 254 62206
## Pattern:Speed 4 0.000022 0.0856 8 62206
## Pattern:as.factor(Channel) 254 0.015139 0.9340 508 62206
## Speed:as.factor(Channel) 254 0.070945 4.5034 508 62206
## Pattern:Speed:as.factor(Channel) 508 0.019778 0.6115 1016 62206
## Residuals 31103
## Pr(>F)
## AgeDays 0.8750
## Pattern 0.9998
## Speed 0.9969
## as.factor(Channel) <2e-16 ***
## Pattern:Speed 0.9996
## Pattern:as.factor(Channel) 0.8534
## Speed:as.factor(Channel) <2e-16 ***
## Pattern:Speed:as.factor(Channel) 1.0000
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
This analysis revealed non-significant effects of AgeDays, Pattern, and Speed, but a significant Channel effect, \(F(254,31103)=15.9913, p<2e-16\), qualified by a significant Speed by Channel interaction, \(F(508,31003)=4.5034, p<2e-16\).
On the basis of these analyses, we drop AgeDays from consideration.
Create mass univariate function.
mass.univariate <- function( ch, df, harm, method ) {
#Select data frame
local.df = df[ df$Channel == ch & df$Harm == harm,]
#Select method and conduct analysis
if (method =="manova.ampl.complex"){
this.stats <- manova( formula = cbind(Sr, Si) ~ Pattern*Speed + Error(iSess), data=local.df ) } else {
this.stats = NA
}
return( this.stats )
}
this.harm = "1F1"
MANOVA on 1F1 with Pattern, Speed, and Pattern x Speed. Analyze and write write output to file.
maov.list = lapply( unique( df.moco.sm$Channel ), mass.univariate, df=df.moco.sm, harm=this.harm, method="manova.ampl.complex")
maov.summ = lapply( maov.list, summary )
Function to extract \(p\) values from MANOVA summary.
# Must customize based base on specific model.
moco.pvals = function( model.list, method="manova.ampl.complex" ) {
ml = unlist( model.list)
# Extract p-values for Pattern, Speed, Pattern*Speed in order
if (method == "manova.ampl.complex") {
new.vals = as.numeric( c( ml['Error: Within.stats21'], ml['Error: Within.stats22'], ml['Error: Within.stats23']) )
}
return( rbind( new.vals ) )
}
Extract \(p\) values from MANOVA summary. Do not adjust for multiple comparisons here as we use our own multiple comparison adjustment procedure. Create data frame for plotting and tweak factor levels and labels.
# Select pvals from MANOVA data
list.pvals = t( sapply(maov.summ, moco.pvals ) )
# Extract p values and create new data frame for plotting. Use unadjusted p values here, as we use our own
adj.Pvals = c(p.adjust(list.pvals[,1], method="none"), p.adjust(list.pvals[,2], method="none"), p.adjust(list.pvals[,3], method="none"))
# Create data frame from pvals for plotting
df.pvals = data.frame( Chan = rep( 1:128, 3),
Cond = rep(c("Pattern", "Speed", "Patt*Spd"), c(128, 128, 128)),
Pvals = as.numeric( adj.Pvals ),
xpos=rep(df.egi$xpos,3),
ypos=rep(df.egi$ypos,3)
)
pvals.cuts = c(-.01,.0001, .0005, .001, .005, .01, .05, 1)
pvals.lbls = c("<.0001","<.0005", "<.001", "<.005", "<.01", "<.05", "ns")
# Create cuts based on p-value levels
Pvals.cuts = cut( df.pvals$Pvals, breaks=pvals.cuts, labels=pvals.lbls)
Pvals.cuts = ordered( Pvals.cuts, levels = rev( pvals.lbls ) )
df.pvals$Pvals.cuts = Pvals.cuts
df.pvals$Cond = ordered( df.pvals$Cond, levels=c("Pattern", "Speed", "Patt*Spd"))
df.pvals$Harm = rep( this.harm, 3*128 )
# Delete scales from plots
yquiet = scale_y_continuous("", breaks=NULL)
xquiet = scale_x_continuous("", breaks=NULL)
#
# # Load topo ears & nose from text file
# img <- readPNG("topoplot.png")
#
# # Add a raster
# topo_ears_nose <- rasterGrob(img, interpolate=TRUE)
pl.title <- paste(this.harm, "Effects by Channel", sep=" ")
pl.theme.topo <- theme( plot.title = element_text(lineheight=.8, face ="bold", vjust=2, size=16),
legend.title=element_text(size=0),
legend.text=element_text(size=12, face="bold"),
legend.position="bottom"
)
# Plot
pl <- ggplot(data=df.pvals, aes(x=xpos, y=ypos, color=Pvals.cuts, size=Pvals.cuts)) +
geom_point() +
facet_grid(facets = . ~ Cond) +
coord_fixed() +
xquiet +
yquiet +
pl.theme.topo
# Add ears, nose
pl + annotation_custom(topo_ears_nose, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf)
Notice right lateral cluster of channels.
chans.below.thresh = df.pvals$Chan[df.pvals$Pvals <= p.thresh ]
df.below.thresh <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == this.harm )
Calculate mean values for each participant, by channel and pattern for the real (Sr) and imaginary (Si) components at the selected harmonic. Then, calculate the group means for Sr and Si and the group standard deviations and root-mean-square (RMS) amplitudes (and standard errors) from the group mean Sr and Si values.
df.below.thresh.pattern <- df.below.thresh %>%
group_by( Channel, Pattern, iSess ) %>%
summarise( sr.sub.mean=mean(Sr),
si.sub.mean=mean(Si),
rms.amp.sub=sqrt(sr.sub.mean^2+si.sub.mean^2)) %>%
group_by( Channel, Pattern ) %>%
summarise( sr.group.mean=mean( sr.sub.mean ),
si.group.mean=mean( si.sub.mean ),
sr.group.sd=sd(sr.sub.mean),
si.group.sd=sd(si.sub.mean),
rms.amp=sqrt( sr.group.mean^2 + si.group.mean^2),
rms.amp.sem = mean( rms.amp.sub )/sqrt(n()),
nsubs=n(),
sr.group.sem=mean(sr.sub.mean)/sqrt(n()),
si.group.sem=mean(si.sub.mean)/sqrt(n())
)
Then plot.
limits <- aes( ymax = rms.amp + rms.amp.sem, ymin = rms.amp )
dodge <- position_dodge( width=0.8 )
pl.theme.bar <- theme(plot.title = element_text(lineheight=.8, face ="bold", vjust=2, size = 20),
panel.background = element_rect(fill=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(fill=NA,color="black", size=.8,
linetype="solid"),
axis.title.x=element_text(vjust=-.6, size=18),
axis.title.y=element_text(face="bold",vjust=1, size=18),
axis.text=element_text(color="black", size=16),
legend.title=element_blank(),
legend.text=element_text(size=16),
legend.position="bottom",
legend.background=element_blank())
pl.pattern <- ggplot( data=df.below.thresh.pattern ) +
aes( x=as.factor(Channel), y=rms.amp, fill=Pattern) +
geom_bar( stat="identity", width=0.8, position=dodge ) +
geom_errorbar( limits, position=dodge, width=0.15 ) +
xlab("Channel") +
ylab(expression(paste("RMS amplitude (", mu, "V)", sep=""))) +
pl.theme.bar +
annotate("rect", xmin = 9.5, xmax = 11.5, ymin = 0, ymax = 1.0, alpha = .1)
pl.pattern
df.below.thresh.speed <- df.below.thresh %>%
group_by( Channel, Speed, iSess ) %>%
summarise( sr.sub.mean=mean(Sr),
si.sub.mean=mean(Si),
rms.amp.sub=sqrt(sr.sub.mean^2+si.sub.mean^2)) %>%
group_by( Channel, Speed ) %>%
summarise( sr.group.mean=mean( sr.sub.mean ),
si.group.mean=mean( si.sub.mean ),
sr.group.sd=sd(sr.sub.mean),
si.group.sd=sd(si.sub.mean),
rms.amp=sqrt( sr.group.mean^2 + si.group.mean^2),
rms.amp.sem = mean( rms.amp.sub )/sqrt(n()),
nsubs=n(),
sr.group.sem=mean(sr.sub.mean)/sqrt(n()),
si.group.sem=mean(si.sub.mean)/sqrt(n())
)
pl.speed <- ggplot( data=df.below.thresh.speed ) +
aes( x=as.factor(Channel), y=rms.amp, fill=Speed) +
geom_bar( stat="identity", width=0.8, position=dodge ) +
geom_errorbar( limits, position=dodge, width=0.15 ) +
xlab("Speed") +
ylab(expression(paste("RMS amplitude (", mu, "V)", sep=""))) +
pl.theme.bar +
annotate("rect", xmin = 9.5, xmax = 11.5, ymin = 0, ymax = 1.0, alpha = .1)
pl.speed
Do for Pattern then Speed. Ignore interaction because there are no effects.
# Summarize by subject first
df.below.thresh.pattern.bysub <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == "1F1" ) %>%
group_by( iSess, Channel, Pattern ) %>%
summarise( sr.mean.bysub=mean(Sr),
si.mean.bysub=mean(Si),
sr.sem=sd(Sr)/sqrt(n()),
si.sem=sd(Si)/sqrt(n())
)
# Then, summarize averages across subjects
df.below.thresh.pattern.summ <- df.below.thresh.pattern.bysub %>%
group_by( Channel, Pattern ) %>%
summarise( sr.mean = mean( sr.mean.bysub ),
si.mean = mean( si.mean.bysub ),
sr.sem = sd(sr.mean.bysub)/sqrt(n()),
si.sem = sd(si.mean.bysub)/sqrt(n())
)
df.below.thresh.speed.bysub <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == this.harm ) %>%
group_by( iSess, Channel, Speed ) %>%
summarise( sr.mean.bysub=mean(Sr),
si.mean.bysub=mean(Si),
sr.sem=sd(Sr)/sqrt(n()),
si.sem=sd(Si)/sqrt(n())
)
df.below.thresh.speed.summ <- df.below.thresh.speed.bysub %>%
group_by( Channel, Speed ) %>%
summarise( sr.mean = mean( sr.mean.bysub ),
si.mean = mean( si.mean.bysub ),
sr.sem = sd(sr.mean.bysub)/sqrt(n()),
si.sem = sd(si.mean.bysub)/sqrt(n())
)
pl.theme.vect <- theme(plot.title = element_text(lineheight=.8, face ="bold", vjust=2, size = 20),
axis.title.x=element_text(vjust=-.6, size=18),
axis.title.y=element_text(face="bold",vjust=1, size=18),
axis.text=element_text(color="black", size=14),
legend.title=element_blank(),
legend.text=element_text(size=16),
legend.background=element_blank(),
legend.position="bottom",
legend.title = element_blank(),
strip.text = element_text( size = 14 ))
df.below.thresh.pattern.summ %>% filter( Channel %in% c(106,107) ) %>%
ggplot() +
aes( x=sr.mean, y=si.mean, color=Pattern ) +
geom_point() +
geom_segment(aes( xend=0, yend=0, x=sr.mean, y=si.mean)) +
geom_pointrange(aes(ymin=si.mean-si.sem, ymax=si.mean+si.sem)) +
geom_errorbarh( aes(xmin=sr.mean-sr.sem, xmax=sr.mean+sr.sem), height=0) +
scale_x_continuous(limits = c(-.6,.6)) +
scale_y_continuous(limits = c(-.6,.6)) +
coord_fixed( ratio=1 ) +
xlab( expression(paste("Signal Real (", mu, "V)", sep=""))) +
ylab( expression(paste("Signal Imaginary (", mu, "V)", sep=""))) +
facet_grid(. ~ Channel) +
pl.theme.vect
pl.theme.vect.all <- theme(axis.title.x=element_text(vjust=-.6, size=18),
axis.title.y=element_text(face="bold",vjust=1, size=18),
axis.text=element_text(color="black", size=8),
legend.title=element_blank(),
legend.text=element_text(size=16),
legend.background=element_blank(),
legend.position="bottom",
legend.title = element_blank(),
strip.text = element_text( size = 14 ))
amp.max = 1.4
df.below.thresh.pattern.summ %>%
ggplot() +
aes( x=sr.mean, y=si.mean, color=Pattern ) +
geom_point() +
geom_segment(aes( xend=0, yend=0, x=sr.mean, y=si.mean)) +
geom_pointrange(aes(ymin=si.mean-si.sem, ymax=si.mean+si.sem)) +
geom_errorbarh( aes(xmin=sr.mean-sr.sem, xmax=sr.mean+sr.sem), height=0) +
coord_fixed( ratio=1 ) +
xlab( expression(paste("Signal Real (", mu, "V)", sep=""))) +
ylab( expression(paste("Signal Imaginary (", mu, "V)", sep=""))) +
facet_wrap(facets= ~ Channel, scales="free") +
pl.theme.vect.all
this.harm="1F2"
MANOVA on 1F2 with Pattern, Speed, and Pattern x Speed. Analyze and write write output to file.
this.harm = "1F2"
maov.list = lapply( unique( df.moco.sm$Channel ), mass.univariate, df=df.moco.sm, harm=this.harm, method="manova.ampl.complex")
maov.summ = lapply( maov.list, summary )
Extract \(p\) values from MANOVA summary. Do not adjust for multiple comparisons here as we use our own multiple comparison adjustment procedure. Create data frame for plotting and tweak factor levels and labels.
# Select pvals from MANOVA data
list.pvals = t( sapply(maov.summ, moco.pvals ) )
# Extract p values and create new data frame for plotting. Use unadjusted p values here, as we use our own
adj.Pvals = c(p.adjust(list.pvals[,1], method="none"), p.adjust(list.pvals[,2], method="none"), p.adjust(list.pvals[,3], method="none"))
# Create data frame from pvals for plotting
df.pvals = data.frame( Chan = rep( 1:128, 3),
Cond = rep(c("Pattern", "Speed", "Patt*Spd"), c(128, 128, 128)),
Pvals = as.numeric( adj.Pvals ),
xpos=rep(df.egi$xpos,3),
ypos=rep(df.egi$ypos,3)
)
pvals.cuts = c(-.01,.0001, .0005, .001, .005, .01, .05, 1)
pvals.lbls = c("<.0001","<.0005", "<.001", "<.005", "<.01", "<.05", "ns")
# Create cuts based on p-value levels
Pvals.cuts = cut( df.pvals$Pvals, breaks=pvals.cuts, labels=pvals.lbls)
Pvals.cuts = ordered( Pvals.cuts, levels = rev( pvals.lbls ) )
df.pvals$Pvals.cuts = Pvals.cuts
df.pvals$Cond = ordered( df.pvals$Cond, levels=c("Pattern", "Speed", "Patt*Spd"))
df.pvals$Harm = rep( this.harm, 3*128 )
pl.title <- paste(this.harm, "Effects by Channel", sep=" ")
# Plot
pl <- ggplot(data=df.pvals, aes(x=xpos, y=ypos, color=Pvals.cuts, size=Pvals.cuts)) +
geom_point() +
facet_grid(facets = . ~ Cond) +
coord_fixed() +
xquiet +
yquiet +
theme( plot.title = element_text(lineheight=.8, face ="bold", vjust=2, size=16),
legend.title=element_text(size=0),
legend.text=element_text(size=12, face="bold"),
legend.position="bottom")
# Add ears, nose
pl + annotation_custom(topo_ears_nose, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf)
Notice midline cluster of channels.
chans.below.thresh = df.pvals$Chan[df.pvals$Pvals <= p.thresh ]
df.below.thresh <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == this.harm )
Calculate mean values for each participant, by channel and pattern for the real (Sr) and imaginary (Si) components at the selected harmonic. Then, calculate the group means for Sr and Si and the group standard deviations and root-mean-square (RMS) amplitudes (and standard errors) from the group mean Sr and Si values.
df.below.thresh.pattern <- df.below.thresh %>%
group_by( Channel, Pattern, iSess ) %>%
summarise( sr.sub.mean=mean(Sr),
si.sub.mean=mean(Si),
rms.amp.sub=sqrt(sr.sub.mean^2+si.sub.mean^2)) %>%
group_by( Channel, Pattern ) %>%
summarise( sr.group.mean=mean( sr.sub.mean ),
si.group.mean=mean( si.sub.mean ),
sr.group.sd=sd(sr.sub.mean),
si.group.sd=sd(si.sub.mean),
rms.amp=sqrt( sr.group.mean^2 + si.group.mean^2),
rms.amp.sem = mean( rms.amp.sub )/sqrt(n()),
nsubs=n(),
sr.group.sem=mean(sr.sub.mean)/sqrt(n()),
si.group.sem=mean(si.sub.mean)/sqrt(n())
)
limits <- aes( ymax = rms.amp + rms.amp.sem, ymin = rms.amp )
dodge <- position_dodge( width=0.8 )
pl.pattern <- ggplot( data=df.below.thresh.pattern ) +
aes( x=as.factor(Channel), y=rms.amp, fill=Pattern) +
geom_bar( stat="identity", width=0.8, position=dodge ) +
geom_errorbar( limits, position=dodge, width=0.15 ) +
xlab("Channel") +
ylab(expression(paste("RMS amplitude (", mu, "V)", sep=""))) +
pl.theme.bar +
annotate("rect", xmin = 5.5, xmax = 6.5, ymin = 0, ymax = 0.5, alpha = .1) + # Chan 62
annotate("rect", xmin = 15-.5, xmax = 15+.5, ymin = 0, ymax = 0.5, alpha = .1) # Chan 75
pl.pattern
df.below.thresh.speed <- df.below.thresh %>%
group_by( Channel, Speed, iSess ) %>%
summarise( sr.sub.mean=mean(Sr),
si.sub.mean=mean(Si),
rms.amp.sub=sqrt(sr.sub.mean^2+si.sub.mean^2)) %>%
group_by( Channel, Speed ) %>%
summarise( sr.group.mean=mean( sr.sub.mean ),
si.group.mean=mean( si.sub.mean ),
sr.group.sd=sd(sr.sub.mean),
si.group.sd=sd(si.sub.mean),
rms.amp=sqrt( sr.group.mean^2 + si.group.mean^2),
rms.amp.sem = mean( rms.amp.sub )/sqrt(n()),
nsubs=n(),
sr.group.sem=mean(sr.sub.mean)/sqrt(n()),
si.group.sem=mean(si.sub.mean)/sqrt(n())
)
pl.speed <- ggplot( data=df.below.thresh.speed ) +
aes( x=as.factor(Channel), y=rms.amp, fill=Speed) +
geom_bar( stat="identity", width=0.8, position=dodge ) +
geom_errorbar( limits, position=dodge, width=0.15 ) +
xlab("Speed") +
ylab(expression(paste("RMS amplitude (", mu, "V)", sep=""))) +
pl.theme.bar +
annotate("rect", xmin = 5.5, xmax = 6.5, ymin = 0, ymax = 0.5, alpha = .1) + # Chan 62
annotate("rect", xmin = 15-.5, xmax = 15+.5, ymin = 0, ymax = 0.5, alpha = .1) # Chan 75
pl.speed
Do for Pattern then Speed. Ignore interaction because there are no effects.
# Summarize by subject first
df.below.thresh.pattern.bysub <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == "1F1" ) %>%
group_by( iSess, Channel, Pattern ) %>%
summarise( sr.mean.bysub=mean(Sr),
si.mean.bysub=mean(Si),
sr.sem=sd(Sr)/sqrt(n()),
si.sem=sd(Si)/sqrt(n())
)
# Then, summarize averages across subjects
df.below.thresh.pattern.summ <- df.below.thresh.pattern.bysub %>%
group_by( Channel, Pattern ) %>%
summarise( sr.mean = mean( sr.mean.bysub ),
si.mean = mean( si.mean.bysub ),
sr.sem = sd(sr.mean.bysub)/sqrt(n()),
si.sem = sd(si.mean.bysub)/sqrt(n())
)
df.below.thresh.speed.bysub <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == this.harm ) %>%
group_by( iSess, Channel, Speed ) %>%
summarise( sr.mean.bysub=mean(Sr),
si.mean.bysub=mean(Si),
sr.sem=sd(Sr)/sqrt(n()),
si.sem=sd(Si)/sqrt(n())
)
df.below.thresh.speed.summ <- df.below.thresh.speed.bysub %>%
group_by( Channel, Speed ) %>%
summarise( sr.mean = mean( sr.mean.bysub ),
si.mean = mean( si.mean.bysub ),
sr.sem = sd(sr.mean.bysub)/sqrt(n()),
si.sem = sd(si.mean.bysub)/sqrt(n())
)
Select channels 62 and 75.
df.below.thresh.speed.summ %>% filter( Channel %in% c(62,75) ) %>%
ggplot() +
aes( x=sr.mean, y=si.mean, color=Speed ) +
geom_point() +
geom_segment(aes( xend=0, yend=0, x=sr.mean, y=si.mean)) +
geom_pointrange(aes(ymin=si.mean-si.sem, ymax=si.mean+si.sem)) +
geom_errorbarh( aes(xmin=sr.mean-sr.sem, xmax=sr.mean+sr.sem), height=0) +
scale_x_continuous(limits = c(-.6,.6)) +
scale_y_continuous(limits = c(-.6,.6)) +
coord_fixed( ratio=1 ) +
xlab( expression(paste("Signal Real (", mu, "V)", sep=""))) +
ylab( expression(paste("Signal Imaginary (", mu, "V)", sep=""))) +
facet_grid(. ~ Channel) +
pl.theme.vect
amp.max = 0.5
df.below.thresh.speed.summ %>%
ggplot() +
aes( x=sr.mean, y=si.mean, color=Speed ) +
geom_point() +
geom_segment(aes( xend=0, yend=0, x=sr.mean, y=si.mean)) +
geom_pointrange(aes(ymin=si.mean-si.sem, ymax=si.mean+si.sem)) +
geom_errorbarh(aes(xmin=sr.mean-sr.sem, xmax=sr.mean+sr.sem), height=0) +
coord_fixed( ratio=1 ) +
xlab( expression(paste("Signal Real (", mu, "V)", sep=""))) +
ylab( expression(paste("Signal Imaginary (", mu, "V)", sep=""))) +
facet_wrap(facets= ~ Channel, scales = "free") +
pl.theme.vect.all
this.harm="2F1"
MANOVA on 2F1 with Pattern, Speed, and Pattern x Speed. Analyze and write write output to file.
maov.list = lapply( unique( df.moco.sm$Channel ), mass.univariate, df=df.moco.sm, harm=this.harm, method="manova.ampl.complex")
maov.summ = lapply( maov.list, summary )
Extract \(p\) values from MANOVA summary. Do not adjust for multiple comparisons here as we use our own multiple comparison adjustment procedure. Create data frame for plotting and tweak factor levels and labels.
# Select pvals from MANOVA data
list.pvals = t( sapply(maov.summ, moco.pvals ) )
# Extract p values and create new data frame for plotting. Use unadjusted p values here, as we use our own
adj.Pvals = c(p.adjust(list.pvals[,1], method="none"), p.adjust(list.pvals[,2], method="none"), p.adjust(list.pvals[,3], method="none"))
# Create data frame from pvals for plotting
df.pvals = data.frame( Chan = rep( 1:128, 3),
Cond = rep(c("Pattern", "Speed", "Patt*Spd"), c(128, 128, 128)),
Pvals = as.numeric( adj.Pvals ),
xpos=rep(df.egi$xpos,3),
ypos=rep(df.egi$ypos,3)
)
pvals.cuts = c(-.01,.0001, .0005, .001, .005, .01, .05, 1)
pvals.lbls = c("<.0001","<.0005", "<.001", "<.005", "<.01", "<.05", "ns")
# Create cuts based on p-value levels
Pvals.cuts = cut( df.pvals$Pvals, breaks=pvals.cuts, labels=pvals.lbls)
Pvals.cuts = ordered( Pvals.cuts, levels = rev( pvals.lbls ) )
df.pvals$Pvals.cuts = Pvals.cuts
df.pvals$Cond = ordered( df.pvals$Cond, levels=c("Pattern", "Speed", "Patt*Spd"))
df.pvals$Harm = rep( this.harm, 3*128 )
pl.title <- paste(this.harm, "Effects by Channel", sep=" ")
# Plot
pl <- ggplot(data=df.pvals, aes(x=xpos, y=ypos, color=Pvals.cuts, size=Pvals.cuts)) +
geom_point() +
facet_grid(facets = . ~ Cond) +
coord_fixed() +
xquiet +
yquiet +
pl.theme.topo
# Add ears, nose
pl + annotation_custom(topo_ears_nose, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf)
Small cluster in right mid-frontal channels for Pattern by Speed interaction, but ps <.001 and p<.01. So, we ignore.
this.harm = "3F1"
MANOVA on 3F1 with Pattern, Speed, and Pattern x Speed.
maov.list = lapply( unique( df.moco.sm$Channel ), mass.univariate, df=df.moco.sm, harm=this.harm, method="manova.ampl.complex")
maov.summ = lapply( maov.list, summary )
Extract \(p\) values from MANOVA summary. Do not adjust for multiple comparisons here as we use our own multiple comparison adjustment procedure. Create data frame for plotting and tweak factor levels and labels.
# Select pvals from MANOVA data
list.pvals = t( sapply(maov.summ, moco.pvals ) )
# Extract p values and create new data frame for plotting. Use unadjusted p values here, as we use our own
adj.Pvals = c(p.adjust(list.pvals[,1], method="none"), p.adjust(list.pvals[,2], method="none"), p.adjust(list.pvals[,3], method="none"))
# Create data frame from pvals for plotting
df.pvals = data.frame( Chan = rep( 1:128, 3),
Cond = rep(c("Pattern", "Speed", "Patt*Spd"), c(128, 128, 128)),
Pvals = as.numeric( adj.Pvals ),
xpos=rep(df.egi$xpos,3),
ypos=rep(df.egi$ypos,3)
)
pvals.cuts = c(-.01,.0001, .0005, .001, .005, .01, .05, 1)
pvals.lbls = c("<.0001","<.0005", "<.001", "<.005", "<.01", "<.05", "ns")
# Create cuts based on p-value levels
Pvals.cuts = cut( df.pvals$Pvals, breaks=pvals.cuts, labels=pvals.lbls)
Pvals.cuts = ordered( Pvals.cuts, levels = rev( pvals.lbls ) )
df.pvals$Pvals.cuts = Pvals.cuts
df.pvals$Cond = ordered( df.pvals$Cond, levels=c("Pattern", "Speed", "Patt*Spd"))
df.pvals$Harm = rep( this.harm, 3*128 )
pl.title <- paste(this.harm, "Effects by Channel", sep=" ")
# Plot
pl <- ggplot(data=df.pvals, aes(x=xpos, y=ypos, color=Pvals.cuts, size=Pvals.cuts)) +
geom_point() +
facet_grid(facets = . ~ Cond) +
coord_fixed() +
xquiet +
yquiet +
pl.theme.topo
# Add ears, nose
pl + annotation_custom(topo_ears_nose, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf)
Interesting midline cluster at p<.0001 and p<.0005 for speed.
chans.below.thresh = df.pvals$Chan[df.pvals$Pvals <= p.thresh ]
df.below.thresh <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == this.harm )
df.below.thresh.speed <- df.below.thresh %>%
group_by( Channel, Speed, iSess ) %>%
summarise( sr.sub.mean=mean(Sr),
si.sub.mean=mean(Si),
rms.amp.sub=sqrt(sr.sub.mean^2+si.sub.mean^2)) %>%
group_by( Channel, Speed ) %>%
summarise( sr.group.mean=mean( sr.sub.mean ),
si.group.mean=mean( si.sub.mean ),
sr.group.sd=sd(sr.sub.mean),
si.group.sd=sd(si.sub.mean),
rms.amp=sqrt( sr.group.mean^2 + si.group.mean^2),
rms.amp.sem = mean( rms.amp.sub )/sqrt(n()),
nsubs=n(),
sr.group.sem=mean(sr.sub.mean)/sqrt(n()),
si.group.sem=mean(si.sub.mean)/sqrt(n())
)
pl.speed <- ggplot( data=df.below.thresh.speed ) +
aes( x=as.factor(Channel), y=rms.amp, fill=Speed) +
geom_bar( stat="identity", width=0.8, position=dodge ) +
geom_errorbar( limits, position=dodge, width=0.15 ) +
xlab("Speed") +
ylab(expression(paste("RMS amplitude (", mu, "V)", sep=""))) +
pl.theme.bar
pl.speed
Do for Pattern then Speed. Ignore interaction because there are no effects.
# Summarize by subject first
df.below.thresh.pattern.bysub <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == "1F1" ) %>%
group_by( iSess, Channel, Pattern ) %>%
summarise( sr.mean.bysub=mean(Sr),
si.mean.bysub=mean(Si),
sr.sem=sd(Sr)/sqrt(n()),
si.sem=sd(Si)/sqrt(n())
)
df.below.thresh.pattern.summ <- df.below.thresh.pattern.bysub %>%
group_by( Channel, Pattern ) %>%
summarise( sr.mean = mean( sr.mean.bysub ),
si.mean = mean( si.mean.bysub ),
sr.sem = sd(sr.mean.bysub)/sqrt(n()),
si.sem = sd(si.mean.bysub)/sqrt(n())
)
df.below.thresh.speed.bysub <- df.moco.sm %>%
filter( Channel %in% chans.below.thresh, Harm == this.harm ) %>%
group_by( iSess, Channel, Speed ) %>%
summarise( sr.mean.bysub=mean(Sr),
si.mean.bysub=mean(Si),
sr.sem=sd(Sr)/sqrt(n()),
si.sem=sd(Si)/sqrt(n())
)
df.below.thresh.speed.summ <- df.below.thresh.speed.bysub %>%
group_by( Channel, Speed ) %>%
summarise( sr.mean = mean( sr.mean.bysub ),
si.mean = mean( si.mean.bysub ),
sr.sem = sd(sr.mean.bysub)/sqrt(n()),
si.sem = sd(si.mean.bysub)/sqrt(n())
)
amp.max = 0.5
df.below.thresh.speed.summ %>%
ggplot() +
aes( x=sr.mean, y=si.mean, color=Speed ) +
geom_point() +
geom_segment(aes( xend=0, yend=0, x=sr.mean, y=si.mean)) +
geom_pointrange(aes(ymin=si.mean-si.sem, ymax=si.mean+si.sem)) +
geom_errorbarh( aes(xmin=sr.mean-sr.sem, xmax=sr.mean+sr.sem), height=0) +
coord_fixed( ratio=1 ) +
xlab( expression(paste("Signal Real (", mu, "V)", sep=""))) +
ylab( expression(paste("Signal Imaginary (", mu, "V)", sep=""))) +
facet_wrap(facets= ~ Channel, scales = "free") +
pl.theme.vect.all
Responses to 4 and 8 deg/s have more similar phase tuning profiles. In fact, phase matters a great deal for channels 84 and 89 where amplitudes look similar.